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Abstract. This paper concerns the propagation of particles through 
a quenched random medium. In the one- and two-dimensional 
models considered, the local dynamics is given by expanding cir- 
cle maps and hyperbolic toral automorphisms, respectively. The 
particle motion in both models is chaotic and found to fluctuate 
about a linear drift. In the proper scaling limit, the cumulative dis- 
tribution function of the fluctuations converges to a Gaussian one 
with system dependent variance while the density function shows 
no convergence to any function. We have verifled our analytical 
results using extreme precision numerical computations. 



1. Introduction 

Variants of a mechanical model now widely known as the Lorentz gas 
have occupied the minds of scientists for more than a century. Initially 
proposed by Lorentz [ ] in 1905 to describe the motion of an electron 
in a metallic crystal, the model consists of fixed, dispersing, scatterers 
in and a free point particle that bounces elastically off the scatterers 
upon collisions. 

If the lattice of scatterers is periodic, the model is also referred to as 
Sinai Billiards after Sinai, who proved [ ] that the system (with d = 2) 
is ergodic if the free path of the particle is bounded. In the latter case 
it was also proved that, in a suitable scaling limit, the motion of the 
particle is Brownian [BuSi]. Sinai's work can be considered the first 
rigorous proof of Boltzmann's Ergodic Hypothesis in a system that 
resembles a real-world physical system. 

Despite its seeming innocence, the Lorentz gas exhibits a great deal 
of complexity. One example is the lack of smoothness of the dynamics 
caused by tangential collisions of the particle with the scatterers. An- 
other one, the presence of recollisions, is a source of serious statistical 
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Figure 1. Schematic diagram illustrating the qualita- 
tive features of our two-dimensional model. The medium 
is composed of two different types of square cells sepa- 
rated by walls which allow particles to pass only in the 
direction of the arrowheads. The zig-zag line shows a 
path of a particle through the medium. 



difficulties that have not been overcome in the study of the aperiodic 
Lorentz gas. For more background, see [la, ibz, UhMa, ChDo] and the 
references therein. 

Our study concerns an idealization of the aperiodic Lorentz gas with 
semipermeable walls illustrated in Figure 1. In each cell, there is a con- 
figuration of scatterers drawn independently from the same probability 
distribution. In our case, the distribution is Bernoulli, so that there 
are two possible configurations to choose from inside each cell. As an 
important aspect, the environment thus obtained is quenched; once the 
scatterer configurations have been randomly chosen, they are frozen for 
good, and the only randomness that remains is in the initial data of 
the particle. Between the cells are semipermeable walls that allow the 
particle to pass through from left to right and from bottom to top, 
as shown by the arrowheads, but not in the opposite directions. This 
model may be thought of as describing the propagation of particles in 
an exotic, anisotropic, medium. 

Notice that, as a significant simplification, there is no recurrence; 
once a particle leaves a cell, it never returns to the same cell again. 
Yet a particle can occupy a single cell for an arbitrarily long time 
before moving on to a neighboring one — albeit a long occupation time 
has a small probability. Moreover, where, when, and in which direction 
the particle exits a cell depends heavily on the scatterer configuration 
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inside the cell, in addition to the position and direction of the particle 
at entry. Inside each cell, the dynamics is chaotic and hyperbolic. 

In our one- and two-dimensional idealizations, the billiard dynamics 
is replaced by discrete dynamical systems acting in each cell. In other 
words, acting on the particle's current position by a map associated 
to the current cell gives its position one time unit later. In dimension 
one the maps associated to the cells are expanding circle maps while in 
dimension two they are hyperbolic toral automorphisms. Such choice 
of maps retains the chaotic and hyperbolic nature of the problem. A 
closely related model has been studied in [\ySt, AyLiSt]. 

Our objective is to understand certain statistical properties of the 
motion. More precisely, we are interested in how the particle dis- 
tribution evolves with time when the initial distribution is uniform 
and supported on one initial cell. We make several analytical proposi- 
tions, which we verify numerically. We show that, on the average, the 
particles follow a linear drift and that, after taking a suitable scaling 
limit, the cumulative distribution of the fluctuations about the mean is 
Gaussian. Moreover, the drift and variance are the same for (almost) all 
environments drawn from the same distribution. Nevertheless, the flne 
structure of the particle distribution is peculiar due to the quenched 
environment. In particular, the density function does not converge to 
that of a normal distribution. In fact, it does not converge at all in the 
aforementioned scaling limit. 

Acknowledgements. We are indebted to Arvind Ayyer and Joel 
Lebowitz for stimulating discussions. Tapio Simula is supported by 
the Japan Society for the Promotion of Science Postdoctoral Fellow- 
ship for Foreign Researchers. Mikko Stenlund is partially supported by 
a fellowship from the Academy of Finland. 

2. One-dimensional model: expanding maps on the circle 

2.1. Introduction. Imagine tiling the nonnegative half line [0, oc) so 
that each interval — or tile — Ik = [k^k + 1) with G N carries a la- 
bel uj{k) that equals either or 1. Such a tiling can be realized by 
flipping a coin for each k and encoding the outcomes in a sequence 
uj = (6j(0), 6^(1), . . . ) G {0, 1}^ called the environment. The coin could 
be balanced but the tosses are independent, with FToh{uj{k) = i) = Pi 
for all k. Below, P^^ will stand for the corresponding probability mea- 
sure on the space of Bernoulli sequences u. In each experiment we 
freeze the environment — meaning that we work with one fixed sequence 
6J at a time. 
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The dynamics in our model is generated by the following definitions. 
Let Ao, Ai G {2, 3, . . . } and define the circle maps Ti{x) = AiX mod 1. 
An experiment comprises iterating many times the map M x §^ ^ 
R X §^ : (v^x) {v + 74^([^])X — x, T^([^])(x)), where [v] is the integer 
part of V and <^([^]) is the corresponding component of u. The actual 
initial condition is (x, x), with x G [0, 1), although in the treatment be- 
low we will also consider initial conditions of the form {V + x^x) with 
y G N arbitrary. Let P denote the Lebesgue measure (i.e., the uni- 
form probability distribution) on the circle §^ and E the corresponding 
expectation. 

The difference A^([^])X — x above can be thought of as the jump 
experienced by a particle. In other words, the evolution of the in- 
coordinate on the line keeps track of the journey of a particle, from the 
point of view of the particle. The x-coordinate is just the projection 
onto the circle (and is needed for technical purposes): if x^) is the 
image of (vq^ Xq) = (1/ + x, x) after n steps, then x^ = Vn mod L 

The model thus describes the deterministic motion of a particle in 
a randomly chosen, but fixed, environment. In probability jargon, the 
particle performs a deterministic walk in a quenched random environ- 
ment. The challenge is that the map determining the position of the 
particle in the nearest future depends on the tile /[^j the particle is in 
through the label cc;(['u]) of the tile. That is, the motion of the particle 
is guided by the a priori chosen environment. 

Example 1. A concrete example is obtained by choosing Aq — 2^Ax — 
3, and Po = Pi = 

In general, 

n-l 

Notice that already in the expression of the initial V and x and the 
environment ou are present in a complicated fashion: 

The complication referred to earlier is obvious; for each x (and cj) we 
have to compute which map the symbol Acc;([^;i]) stands for. Each Vn 
{n > 2) is generically a piecewise affine function of x and the number 
of discontinuities grows exponentially with n. 

We anticipate that, for large values of n, behaves statistically (in 
the weak sense) like 

M{nD,na'^), (1) 
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where D is a deterministic number called the drift and jV{nD^na'^) 
stands for a real-valued, normally distributed, random variable with 
mean nD and variance na^. In principle, D and could depend on the 
environment cj, but remarkably it turns out that they do not, as long 
as the environment is typical By typicality we mean that u belongs 
to a set whose P^q -probability is one and whose elements enjoy good 
statistical properties such as the convergence of ^#{k < n \ uj{k) = 0} 
to the limit po. 

It is reasonable to expect that the limit 

exists and has the same value for almost all x ^. Thus, we are led to 
conclude that 

= lim = lim -EK(x)) . 



The final equality follows from the bounded convergence theorem. 
Let us consider the (asymptotically) centered random variable 



which measures the fluctuations of relative to the linear drift. Pro- 
vided (1) is true, is approximately Gaussian with variance na^ 
More precisely, we would like to know if -^X^ converges in distribu- 
tion to A/'(0, cr^). By deflnition, this means that, for any flxed ?/ G M, 

lim p(^Xn <y] = r e"''/'"' ds. 

n^oo J ^/27^aJ-oo 



2.2. Markov partition. We next reduce the deterministic walk in a 
random environment to a random walk in a random (still quenched) 
environment which is easier to treat. This can be done using a Markov 
property of the tiling that allows us, in the statistical sense, to ignore 
the exact position of the particle and only keep track of the tile it is 
occupying. 

Let [ • ] denote the integer part of a number. If we deflne 

K = [Vfi] and Xn = Vn- K], (2) 

-'^This cannot hold for all x. For instance, if x = then 1 = Vk = Vk-\-i = . . . 

and the process stops. 
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then the earher dynamics with the initial condition ('^0,^0) = (Vo + 
Xq^Xq) is equivalent to 

K+l = K + [A^{Vn)Xn\ 

Recall our convention [0, oc) = IJ^lo^^' where Ik = [k^k + 1) is 
called a tile. Suppose now that ^ h- This is equivalent to ^4 = fc. 
As before, we are interested in the probability distribution of but 
this time only at the level of tiles. More precisely, we wish to know 
the probability distribution of Vn- This is the probability vector p^^^ = 
{p^^\ p^i \ . . . ) where the numbers 

pi") = P(K = k) 

are such that Yl'k={)Pk^ = ^^ 

We now consider the dynamical system being initialized with the 
condition ('Uo,Xo) = (Vq + Xo,Xo), where Vq G N and G [0,1) are 
random variables; is uniformly distributed and independent of Vq, 
but for the moment we do not specify the distribution pf^ = P(Vo = k) 
which is the initial distribution of the chain Vn- 

Since each Ail^ is exactly the union of a few of the intervals 7/^/^, the 
collection {1^} is a simultaneous Markov partition for the two maps. 
We then obtain the Markov property 

P{Vn = K I Vn-1 = . . . , H = /Cq) = P(K = k^ \ K-1 = k^-l) 

for arbitrary histories. Thus, ^4 is a time-homogeneous Markov chain 
on the countably infinite state space N with the transition probabilities 



PK+i G IkM I ^ h) = if / G {0, 1, . . . , A^^k) - 1}, 

otherwise. 



Notice that the above holds for any environment, cj, but the resulting 
Markov chain does depend on the choice of u. 

Defining the transition matrix F = {'^k^k')kk'^ P^^^ ^ p^'^~^^T. Thus, 

for an arbitrary initial distribution. In particular, 

p(°) = (1,0,0,...) 

corresponds to initializing the dynamical system at {v^^x^)) = 
with < X < 1 arbitrary. 



Ai maps A: + ^, A: + affinely onto [A: + /, A: + / + 1). 
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In principle, (4) provides us with complete statistical understanding 
of the dynamics. For instance, the drift can be expressed as 

1 1 1 ^ 

D = lim - E{vn) = lim - E(K) = lim -V 

In practice, calculating for large values of n is difficult. 

2.3. Drift and variance. For each (i, j) G {0, 1}^ the transition prob- 
ability at time n from a tile labeled i to a tile labeled j is 

aij{n) = P(6j(K+i) = j I ^(K) = i)- 

The analysis of this quantity is subtle, because it depends on the tiling. 
For instance, if = (0, 0, . . . ), then P(6j(K+i) = | uj{Vn) = 0) = 1. 
The conditional probability P x Fpq{uj{Vi) = j \ ^(Vo) = i) equals 

because the elements uj{k) of the tiling are independent. Here pi is the 
Bernoulli probability of getting an i in the tiling. We think of a^j as 
an effective transition probability which only depends on the statistical 
properties of the tiling. 

We expect the actual transition probability aij{n) to converge to the 
effective value a^j with increasing time, 

lim aij{n) = a*-. 

This is so because, as n increases, the position of the particle at time n 
depends on the tiling on an increasing subinterval of [0, oc) and should 
therefore reflect increasingly the statistics of the tiling instead of its 
local details. 
Moreover, 

lim(aT=f^ (5) 

for a p G (0,1) that can be found by diagonalizing a* or by solving the 
equilibrium equation (p, 1 — p)^* = (p, 1 — p): 

l-Pii^-Poi^ AoAi - piAi - poAo' 

For instance, in the case of Example 1 we obtain p = |. 
Notice that, for any probability vector (g, 1 — q), 

{q,l-q) lim = -p). 
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The probability vector 

k 

{q,l-q)Y[a{n) = {q{k),l-q{k)) 
will converge to some (g*, 1 — g*), because a{n) ^ a*. In fact, 

2N 2N 

lim (g, 1 - g) TT a{n) = lim {q{N), 1 - q{N)) TT a{n) 

n=0 n=N+l 

= iq\l-q*)]im ia*f = ipA-p), 

as ^ oc. 

We interpret the result above so that P{uj{Vn) = 0) ^ p and 
P{uu{Vn) = 1) 1 — p as n ^ oc. That is, along a given (typi- 
cal) trajectory, the fraction of time the particle spends in a tile labeled 
is p: 

Um*<^<"l-"^^'=°>=p. (6) 

Notice that p does not depend on the (typical) tiling. 

2.3.1. Dnft Define the jumps = Vi - Vi^i {i > 1). Then K = 
Yl^=i^i' We also denote ^^^"^ a random variable that takes values in 
{0, . . . , Aj — 1} with uniform distribution. For a (typical) tiling, 

D = lim = lim =pE(^(°)) + {1 - p)E{^^^^) 

Ao-l .Ai-l pAo + {l-p)Ai-l 

=p^ + (l-p)^- = ^ • 

We also claim that this equals lim^^oo ^ for almost all x. In the case 
of Example 1 , D = | . 

2.3.2. Variance. The variance is 



= lim Varf — ^ ] = lim ^ — 0^ 



n 



Let us assume < Ai and study the process Wn = ^^^i Ci having 
the i.i.d. increments Q whose distribution is Prob(Ci = k) = ^ + 
if < < ^0 and Prob(Ci = k) = ^ if Aq < k < Ai. The increments 
have been chosen so that Wn mimics Vn as closely as possible and has 
the same variance. For instance, staying in the same tile {(i = 0) has 
probability -j- + where p is the probability of being in a tile labeled 



s the prol 

term accounts similarly for the case of label 1. Then Mean(W^) 



and ^ is the probability of staying in that tile, while the second 
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nMean((^i) nD. Setting K{m) Yll^=^ + + 6^' 

variance of Wn is 

Var(W^„) = n Var(Ci) = n (^^(A) + ^^(^i) " ^ 

Using the values p = | and D = ^ obtained from Example 1, the 
formula above gives ^ Var(W^) = ||. 

2.4. Sensitivity on the initial condition. Let us next consider the 
Lyapunov exponent 

A = lim — m — — = lim — > m 

k=l ^ ^ 

which measures the exponential rate at which two nearby initial points 
drift apart under the dynamics. Above, the chain rule has been used. 
Recall the notation introduced in (2) and that 'Uq = ^0 = ^- As '^^^^^ = 



With probability zero Vk-i is an integer, in which case Vk-i = Vk-i-> 
Xk-i = 0, and the process stops. We assume that Vk-i is not an integer. 

Then '-^^ = 0, = 1, and ^ = A^^y,_,), such that, by (6), 



1 

A lim - y^lnA^fVk_i) =plnAo + (1 -p)lnAi 

and is positive. Roughly speaking, the distance between two very 
nearby trajectories thus grows hke e^" = (A^Al"^) , which is tan- 
tamount to chaos. 

2.5. Numerical study. In order to study the model introduced above 
numerically we first create the random tiling (or environment) uj of 
length 2n+l where n is the number of jumps to be performed in a single 
trajectory. This guarantees that every possible path fits inside the tiling 
although some computational eflFort could be saved by choosing the 
number of tiles closer to [nD]. Notice also that while in principle new 
tiles could be added dynamically to the end of the tiling as required, it 
is computationally far more efficient to construct the tiling as a static 
entity in the beginning of the computation. 

In practice, the tiling is generated by producing a vector of pseudo- 
random numbers distributed uniformly on the interval (0, 1) using the 
Mersenne Twister algorithm. The label of the tile cj(/c) is then obtained 
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n 



Figure 2. Position of a particle Vn divided by the num- 
ber of jumps n taken for a single typical trajectory. 
The horizontal line denotes the exact value for the drift 
D = 5/7. The inset shows a blow-up of a part of the 
main figure. 

by rounding the number on each tile k to the nearest integer. The 
computational tiling u is finalized by the operation u{k) = Lju{k){Ai — 
Aq) + Aq yielding a vector whose each element is either 2 or 3. 

Each ensemble member (particle trajectory) is initialized by generat- 
ing a pseudo-random number to determine the starting point Xq G (0, 1) 
of the trajectory. The subsequent particle positions are determined by 
the underlying tiling. The jumping process could be performed deter- 
ministically by keeping track of the exact position of the particle. 
However, the Markov property of the process provides us a superior way 
of obtaining the desired statistics stochastically. In this algorithm, be- 
fore every jump, we sample a new pseudo-random number d from the 
interval (0,cj(/c)) depending on the current tile k. Then a jump to the 
tile k + [d] is made and the whole procedure is repeated n times to 
produce a single trajectory. 

Figure 2 shows a typical trajectory of n = 10^ jumps obtained using 
the above prescription. The position of the particle Vn divided by the 
number of jumps n taken is clearly seen to saturate to the analytical 
value for the drift D = 5/7^ plotted as a straight line in the figure. The 
inset shows the late-time evolution of the drift of the particle. 

Figure 3 displays the computed probability density function for the 
random variable (Vn—nD) / ^Jn obtained using 10'' trajectories of length 
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{Vn-nD)/^ 



Figure 3. Probability density of the random variable 
{Vn — nD)/^/n. The solid curve is the Gaussian with 
zero mean and variance = 24/49. The insets in the 
left- and right-hand sides show the probability densities 
for the subsets of trajectories ending to a tile labeled by 
and 1, respectively. The horizontal lines at levels 0.4, 
0.8, 0.3 and 0.6 in the insets are plotted to guide the eye. 



n = 10^. The solid curve is a normalized Gaussian function with zero 
mean and variance of = 24/49. Each point in the figure indicates 
the relative frequency that a trajectory ends in the corresponding tile 
after n steps. In the left- and right-hand side insets only trajectories 
ending in tiles labeled by and 1, respectively, are considered. If the 
graph in the right-hand side inset is vertically stretched by the factor 
it becomes practically overlapping with the one in the left-hand 
side inset. In the figure one can discern several Gaussian shapes, all 
of which are very well approximated by the analytical Gaussian after 
normalization with a suitable constant. 

Figure 4 shows two cumulative distribution functions, obtained by 
integrating the numerical and analytical probability densities shown in 
Figure 3. They match to a great accuracy and we are lead to believe 
that the random variable X^/ \/n = {vn — nD) / ^Jn is, indeed, normally 
distributed with zero mean and variance . We have also analyzed 
the characteristic function which leads to the same conclusion. 
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0^— =^ ^ ^ ^ 1 

-2-1012 



Figure 4. Cumulative distribution function corre- 
sponding to the data shown in Figure 3. The two curves 
shown (computational and analytical) are overlapping. 

To conclude, our numerical data strongly supports the theoretical 
analysis presented earlier. Within the numerical accuracy, the distri- 
bution is Gaussian with the drift and variance predicted by our ana- 
lytical calculations. We have performed these numerical experiments 
using different (fixed) tilings and filling probabilities, p^, and have al- 
ways arrived at the same conclusion. 

3. TW^O-DIMENSIONAL MODEL: HYPERBOLIC TORAL 
AUTOMORPHISMS 

3.1. Introduction. We begin by tiling the first quadrant of the plane 
by unit squares, attaching the label or 1 to each tile. That is, corre- 
sponding to each vector = (/ci, G the tile [/ci, fci + l) x [/c2, ^2 + 1) 
carries a label 6j(/c) G {0,1}. The fixed tiling bj = (6j(/c))/eGN2 is our 
environment and P^q stands for the Bernoulli probability measure on 
the space of such tilings. 

The process takes place on the plane and each Ai {i = 0, 1) 
is a matrix with positive integer entries and determinant 1. Such a 
matrix is hyperbolic, with two eigenvalues, A > 1 and A~^, and the 
eigenvector corresponding to A points into the first quadrant. The 
formula TiX = Aix mod 1 defines a hyperbolic toral automorphism. 
A precise description of the dynamics is given by the map x ^ 
X : {v^x) 1-^ ('^ + — x, T^([^])(x)), where [v] is the integer 
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part of V. The initial condition is (x, x), with x G [0, 1)^. Let P denote 
the Lebesgue measure (i.e., the uniform probabihty distribution) on 
the torus and E the corresponding expectation. 

Example 2. A concrete example is obtained by choosing Aq = (i j); 
Ai = (ID, andpQ =pi = |. 

We claim that the limit 



D = e( lim = lim -EK(x)), 

called the drift, exists and that the (asymptotically) centered random 
vector 

which measures the fluctuations of Vn relative to the linear drift, is 
approximately Gaussian with covariance matrix na^. More precisely, 
-^Zn converges in distribution to A/'(0,a^), where is given by 

lim^^ooCov(^;^Z^, denoting = (-oc,Zi] x (-00,^2] for 

any flxed z = (zi, Z2) G M^, 

lim P (^X^ < zi, -^Y^ < Z2) = ^L= / e-^-(-^)-^^ d's. 



In contrast with the one-dimensional case, the tiling is not a Markov 
partition for the maps, which considerably complicates the analysis of 
the model. 

3.2. Drift. Let us continue to denote Vn = [vn]- We conjecture that 
the drift vector D = ( ) is given by 

D = pDo + {l-p)Di, 

where Di = Ei{AiX — x) = {Ai — 1) ^ | ^ equals the average jump 

under the action of the matrix Ai and p is as in (6). The value of 
p is obtained, as above (5), from an eflFective transition matrix a*. 
Its general element a*^- is the conditional probability P x ¥pq{uj{Vi) = 
j I u{Vo) =i) = P X ¥p^{uj{[Aix]) = j 1 6j(0, 0) = i)^the probability of 
jumping to a tile labeled j when the initial tile is labeled i and when 
the choice of the tiling is being averaged out. 

In practice, a*^- is computed as follows. We assume that the initial tile 
is labeled i.e., 6j(0, 0) = i. The image of the unit square under A^ is a 
parallelogram of area one that overlaps with various tiles. The area of 
intersection of the parallelogram with a tile represents the probability 
of jumping to that tile, a*^- can then be computed recalling that each 
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tile is labeled with probability po independently of the others. In 
the case of Example 2, we obtain Dq = (^l^ ^ Di = (^^^ ^ and a* = 

(1^31 31 \ / 5 3 \ / \ 

"It 2 1^1 1 = 1 i , which results in p = § and = f| . 
6 2 6^6 2 / V 12 12 / V 19 / 

3.3. Numerical study. As mentioned earlier, the Markov property 
deployed in the numerical study of the one-dimensional problem where 
we used a stochastic jumping algorithm does not, unfortunately, apply 
in the two-dimensional case. Instead, we are forced to compute the 
particle trajectories fully deterministically which renders the numeri- 
cal problem difficult. Due to the chaotic nature of the process, the 
position Vn of the particle must now be represented with an accuracy 
to approximately 2n decimal places in order to keep the accumulation 
of the numerical rounding errors bounded. This must be done using a 
software implementation since the double precision float native to the 
hardware only contains 15 decimal places. 

We first create the tiling u as in the one-dimensional case with the 
exception that it is now a two-dimensional object. We then choose 
randomly the initial position of the particle within the unit square. 
The label 6j(0,0) of the initial tile is then read and the new position 
of the particle is computed by applying the corresponding map T^{[v])- 
This jumping procedure is repeated n times. The time to compute 
a single trajectory increases dramatically as the path length n is in- 
creased due to the corresponding increase in the required accuracy of 
the representation of the position of the particle. 

Figure 5 shows the convergence of the drifts Di and that of the 
covariance matrix elements afj. The values in decending order at n = 
10^ are ^1,^2, cr^i, (^12 = ^21 5 0-22. The straight lines indicate the 
analytical values for the drift components. Each data point is composed 
using 10^ trajectories. 

In Figure 6 (a)-(c) we have plotted the particle positions in the plane 
after n = 1,2,3, and 2000 jumps, respectively. The trajectories were 
initiated randomly from the unit square. The straight diagonal lines 
indicate the direction of the drift and the cross in frame (d) denotes 
the directions of the eigenvectors of the covariance matrix. Since the 
initial tile in our environment had the label 6j(0,0) = 1, the frame 
(a) simply shows how Ai maps the unit square. The subsequent jumps 
shred the distribution, as illustrated by the frames (b) and (c), because 
particles in different tiles undergo different transformations. Figure 7 
shows a contour plot of the particle distribution after n = 100 jumps 
and reveals prominent stripes, due to the shredding, which are roughly 
aligned with the direction of the drift vector. 
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Figure 5. Values for the x and y components of the 
drift and the covariance matrix elements af^^ cr^2 = ^215 
and (7 22 as a function of n, respectively, in decending 
order at n = 1000. The straight lines indicate the ana- 
lytical values of drift components di and di. 




1 2 3 0123456 

X X 




X X 



Figure 6. End points of 10^ trajectories in the plane 
after n = 1 (a), n = 2 (b), n = 3 (c), and n = 2000 (d) 
jumps. The straight diagonal lines trace the the drift 
vector and the cross in frame (d) shows the eigendirec- 
tions of the covariance matrix. Each frame comprises 
10000 data points. 
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Figure 7. Contour plot of the particle distribution in 
the plane after n = 100 jumps. The straight line shows 
the direction of the drift vector. 

-3 

X 10 



1^ 




Figure 8. Probability density function for the random 
variable Z^/ \/n. The spikes are an inherent feature of 
the distribution and the density does not converge to any 
function. Embedded is the analytical Gaussian function. 

Figure 8 shows the probability density of -^Z^ obtained after n = 
2000 jumps. Embedded is also a two-dimensional Gaussian probability 
density which has the same covariance matrix as the numerical data. 
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Figure 9. Cumulative distribution function corre- 
sponding to the data shown in Figure 8. 

Despite of the fact that the density function itself does not converge 
to any function, the corresponding cumulative distribution function 
shown in Figure 9 is smooth and matches that of the corresponding 
Gaussian distribution. The maximum absolute difference between the 
numerical and analytical functions is 0.017, most of which is due to the 
highest peak in Figure 8. 



4. Conclusions 

We have investigated the statistical properties of a deterministic walk 
in a quenched one-dimensional random environment of expanding cir- 
cle maps and have analytically found the drift and variance for the 
resulting Gaussian probability distribution. Using numerical experi- 
ments we have been able to verify our analytical predictions. We have 
further studied a two-dimensional model similar to the one-dimensional 
system where hyperbolic toral automorphisms take the place of the cir- 
cle maps. Again the probability distribution turns out to be Gaussian 
with certain linear drift and variance. The key feature and compli- 
cating factor in both the one- and two-dimensional cases is the fixed 
random environment. A direct consequence of this is that, even after 
the proper scaling, the probability density does not converge to any 
function — a result which persists both in our one- and two-dimensional 
models. The implementation of recurrence to our model will be left for 
future work. 
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